import numpy as np
import pandas as pd
from sympy import symbols, solve
from nodepy import rk
# Visualisation libraries
## Text
from colorama import Fore, Back, Style
from IPython.display import Image, display, Markdown, Latex, clear_output
## plotly
from plotly.offline import init_notebook_mode, iplot
import plotly.graph_objs as go
import plotly.offline as py
from plotly.validators.scatter.marker import SymbolValidator
from plotly.subplots import make_subplots
# Graphics in retina format
# %config InlineBackend.figure_format = 'retina'
import warnings
warnings.filterwarnings('ignore')
In this section, we go through several methods that can approximate the solution of a well-posed initial-value problem (IVP) on grid points. Consider the following IVP:
\begin{align} \frac{dy}{dt} &= f(t,y),& a \leq t \leq b,~y(a) = y_{0}. \end{align}Here it is assumed that the grid points are equally distributed throughout the interval $[a, b]$. Let $N$ be a positive integer, we have, the following step size
\begin{align} h = \Delta t = \frac{b-a}{N} \end{align}and each grid point can be identified as \begin{align} t_{j} &= a +jh,& j = 0,1,2,\ldots,N. \end{align}
Moreover, let $y_{j}$ denote approximation $y$ at grid points for $j = 1,2,3,\ldots,N$, we can use several methods.
In this article, we only discuss explicit methods Runge–Kutta methods. Explicit Runge-Kutta methods take the form \begin{align} y_{n+1} &= y_{n}+h\sum _{i=1}^{s}b_{i}k_{i},\\ k_{i} &= f\left(t_{n}+c_{i}h,y_{n}+h\sum _{j=1}^{i-1}a_{ij}k_{j}\right). \end{align}
This process is carried forward iteratively until point $x_{n-1}$ is reached. To represent the coefficients of this method, the Butcher tableau is often used which puts the coefficients of the method in a table as follows.
\begin{align} \begin{array}{c|cccc} c_1 & a_{11} & a_{12}& \dots & a_{1s}\\ c_2 & a_{21} & a_{22}& \dots & a_{2s}\\ \vdots & \vdots & \vdots& \ddots& \vdots\\ c_s & a_{s1} & a_{s2}& \dots & a_{ss} \\ \hline & b_1 & b_2 & \dots & b_s\\ \end{array} \end{align}Notes:
Most of the previous explicit methods can be obtained using the Butcher tableau.
This method is defined as follows, \begin{align} y_{n+1} &= y_{n}+\frac{h}{6}\left( k_{1} + 4k_{2} + k_{3} \right),\quad \text{for }n = 0, 1, 2, 3, \ldots, \end{align} with \begin{align} \begin{cases} k_{1} & = f\left(t_{n},~y_{n} \right),\\ k_{2} & = f\left(t_{n}+\frac{h}{2},~y_{n}+\frac{h}{2}k_{1} \right),\\ k_{3} & = f\left(t_{n}+h,~y_{n}+ h(2k_{2} -k_{1}) \right) \end{cases} \end{align} or using the Butcher tableau \begin{align}\begin{array}{c|ccc}0&0&0&0\\1/2&1/2&0&0\\1&-1&2&0\\\hline &1/6&2/3&1/6\\\end{array}\end{align}
Furthermore, we can prepare a Python code using the above algorithm.
def Runge_Kutta_3rd(f, y0, a, b, h= False, N=False):
'''
Inputs:
f: the ODE y'=f(t,y(x))
y0: the initial value
x_range: interval [a,b]
N: Number of points
h: step size h
'''
if N:
h = (b-a)/(N)
if h:
N = int((b-a)/h)
t = np.linspace(a, b, N+1)
y = np.zeros(t.shape, dtype=float)
y[0] = y0
for i in range(N):
k1 = f(t[i], y[i])
k2 = f(t[i]+ h/2, y[i]+ h*k1/2)
k3 = f(t[i]+ h, y[i] + h*(- k1 + 2*k2))
y[i+1] = y[i] + h*(k1 + 4*k2 + k3)/6
del k1, k2, k3
Table = pd.DataFrame({'t':t, 'y':y})
return Table
Example: Conisder the initial value problem $\begin{cases}y'+2ty=te^{-t^2},\\ y(0) = 1 \end{cases},\quad 0 \leq x \leq 1$ with exact solution $y\left(x \right) = \left(1 + \dfrac{t^{2}}{2}\right) e^{- t^{2}}$. We use the the above method for solving this IVP.
# f(t, y(t)):
f = lambda t, y: t*np.exp(-t**(2))-2*t*y
(a, b) = (0, 1)
# the eact solution y(t)
y_exact = lambda t: (1+(t**2)/2)*np.exp(-t**2)
#
y0 = 1
For convenience, we can define the following function which can be used for various methods.
def Figs1(Solver, N = 10, f = f, y0 = y0, a = a, b = b, yLim = [0.4, 1.2], size = [500, 980], Title = False):
Table = Solver(f = f, y0 = y0, a = a, b = b, N = N)
Table['Exact'] = y_exact(Table['t'])
Table['Error'] = np.abs(Table['Exact'] - Table['y'])
display(Table[1:].style.set_properties(subset=['Error'], **{'background-color': 'Lavender', 'color': 'Navy',
'border-color': 'DarkGreen'}).set_precision(4).format({'Error': "{:.4e}"}))
fig = make_subplots(rows=1, cols=2, specs=[[{"type": "xy"}, {"type": "xy"}]], horizontal_spacing = .12)
# Left
fig.add_trace(go.Scatter(x = Table['t'], y = Table['y'], mode='lines', name='Approximation',
line=dict(color='royalblue', width=2)), row=1, col=1)
fig.add_trace(go.Scatter(x = Table['t'], y = Table['Exact'], mode='lines', name='Exact',
line=dict(color='ForestGreen', width=2)), row=1, col=1)
# Right
fig.add_trace(go.Scatter(x = np.arange(Table.shape[0]), y = Table['Error'],mode='markers', name='Error',
marker=dict(color='OrangeRed', size=5, line=dict(color='Black', width=1))), row=1, col=2)
#
fig.update_xaxes(showline=True, linewidth=1, linecolor='Lightgray', mirror=True, zeroline=True, zerolinewidth=1,
zerolinecolor='Black', showgrid=True, gridwidth=1, gridcolor='Lightgray')
fig.update_yaxes(showline=True, linewidth=1, linecolor='Lightgray', mirror=True, zeroline=False, zerolinewidth=1,
zerolinecolor='Black', showgrid=True, gridwidth=1, gridcolor='Lightgray')
# Left
fig.update_xaxes(title_text = '$t$', range= [Table['t'].min(),Table['t'].max()], row=1, col=1)
fig.update_yaxes(title_text = '$y$', range= yLim, row=1, col=1)
# Right
fig.update_xaxes(title_text = 'N', range= [0, Table.index.max()], row=1, col=2)
fig.update_yaxes(title_text = 'Error', type="log", row=1, col=2)
fig.update_layout(legend_orientation='v', plot_bgcolor= 'white', height= size[0], width= size[1])
fig.update_layout(legend=dict(x=1.02, y=.98, font=dict(size=12, color="black"), bordercolor="Black", borderwidth=1))
if Title:
fig.update_layout(title={'text': '<b>' + Title + '<b>','x':0.46, 'y':0.88, 'xanchor': 'center', 'yanchor': 'top'})
fig.show()
def Figs2(Solver, h = [2**(-i) for i in range(3, 12)],
f = f, y0 = y0, a = a, b = b,
yLim = [2.9, 3.1], size = [550, 980], Title = False):
Cols = ['h', 'N', 'Eh']
Table = pd.DataFrame(np.zeros([len(h), len(Cols)], dtype = float), columns=Cols)
Table['h'] = h
Table['N'] = ((b-a)/Table['h']).astype(int)
for n in range(Table.shape[0]):
TB = Solver(f = f, y0 = y0, a = a, b = b, h = Table['h'][n])
Table['Eh'][n] = np.max(np.abs(y_exact(TB['t'])[1:] - TB['y'][1:]))
fig = make_subplots(rows=1, cols=2, horizontal_spacing = 0.02, column_widths=[0.6, 0.4],
specs=[[{"type": "scatter"},{"type": "table"}]])
s = np.log2(Table['Eh'].values[:-1]/Table['Eh'].values[1:])
fig.add_trace(go.Scatter(x = np.arange(1, len(s)+1), y = s, mode='markers',
marker_symbol='circle',
marker_line_color="midnightblue", marker_color="lightskyblue",
marker_line_width=2, marker_size=10), row=1, col=1)
fig.update_xaxes(showline=True, linewidth=1, linecolor='Lightgray', mirror=True,
dtick=1,
# zeroline=True, zerolinewidth=1, zerolinecolor='Black',
showgrid=True, gridwidth=1, gridcolor='Lightgray', title_text = '$h_i$',
range= [0.8, len(s)+.2], row=1, col=1)
fig.update_yaxes(showline=True, linewidth=1, linecolor='Lightgray', mirror=True,
zeroline=True, zerolinewidth=1, zerolinecolor='Black',
showgrid=True, gridwidth=1, gridcolor='Lightgray',
title_text = '$\\log_{2} \\left( \\frac{E_{h_{i}}}{E_{h_{i-1}} }\\right)$', range= yLim, row=1, col=1)
fig.update_layout(plot_bgcolor= 'white', height= size[0] , width= size[1], showlegend=False)
T = Table.copy()
T[['h','Eh']] = T[['h','Eh']].applymap(lambda x: '%.4e' % x)
Temp = []
for i in T.columns:
Temp.append(T.loc[:,i].values)
fig.add_trace(go.Table(header=dict(values = list(Table.columns), line_color='darkslategray',
fill_color='Navy', align=['center','center'],
font=dict(color='white', size=12), height=25), columnwidth = [0.4, 0.2, 0.4],
cells=dict(values=Temp, line_color='darkslategray', fill=dict(color= ['white', 'white', 'Lavender']),
align=['center', 'center'], font_size=12, height=20)), 1, 2)
del Temp, T
if Title:
fig.update_layout(title={'text': '<b>' + 'Order of accuracy (%s)' % Title + '<b>', 'x':0.5,
'y':0.90, 'xanchor': 'center', 'yanchor': 'top'})
else:
fig.update_layout(title={'text': '<b>' + 'Order of accuracy' + '<b>', 'x':0.5,
'y':0.90, 'xanchor': 'center', 'yanchor': 'top'})
fig.show()
def Header(Text, L = 120, C1 = Back.BLUE, C2 = Fore.BLUE):
print(C1 + Fore.WHITE + Style.NORMAL + Text + Style.RESET_ALL + ' ' + C2 +
Style.NORMAL + (L- len(Text) - 1)*'=' + Style.RESET_ALL)
def Line(L=120, C = Fore.BLUE): print(C + Style.NORMAL + L*'=' + Style.RESET_ALL)
Header(Text = 'Numerical Solution using the third-order Runge–Kutta method and N = 10', C1 = Back.GREEN, C2 = Fore.GREEN)
Figs1(Solver = Runge_Kutta_3rd, Title = 'Third-order Runge–Kutta method')
Header(Text = 'Order of Accuracy')
Figs2(Solver = Runge_Kutta_3rd, h = [2**(-i) for i in range(4, 13)], yLim = [2.9, 3.1],
Title = 'Third-order Runge–Kutta method')
Line()
Numerical Solution using the third-order Runge–Kutta method and N = 10 =================================================
| t | y | Exact | Error | |
|---|---|---|---|---|
| 1 | 0.1000 | 0.9950 | 0.9950 | 8.4271e-06 |
| 2 | 0.2000 | 0.9800 | 0.9800 | 1.7039e-05 |
| 3 | 0.3000 | 0.9551 | 0.9551 | 2.5794e-05 |
| 4 | 0.4000 | 0.9203 | 0.9203 | 3.4545e-05 |
| 5 | 0.5000 | 0.8762 | 0.8762 | 4.2967e-05 |
| 6 | 0.6000 | 0.8233 | 0.8233 | 5.0491e-05 |
| 7 | 0.7000 | 0.7628 | 0.7627 | 5.6286e-05 |
| 8 | 0.8000 | 0.6961 | 0.6960 | 5.9311e-05 |
| 9 | 0.9000 | 0.6251 | 0.6250 | 5.8454e-05 |
| 10 | 1.0000 | 0.5519 | 0.5518 | 5.2752e-05 |
Order of Accuracy ======================================================================================================
========================================================================================================================
This method is defined as follows, \begin{align} y_{n+1} &= y_{n}+\frac{h}{4}\left( k_{1} + 3k_{3} \right) \end{align} with \begin{align} \begin{cases} k_{1} & = f\left(t_{n},~y_{n} \right),\\ k_{2} & = f\left(t_{n}+\frac{h}{3},~y_{n}+\frac{h}{3}k_{1} \right),\\ k_{3} & = f\left(t_{n}+\frac{2}{3}h,~y_{n}+ \frac{2h}{3}k_{2}\right) \end{cases} \end{align} or using the Butcher tableau \begin{align} \begin{array}{c|ccc}0&0&0&0\\1/3&1/3&0&0\\2/3&0&2/3&0\\\hline &1/4&0&3/4\\\end{array} \end{align}
The Butcher tableau can be also generated using nodepy python package. For example, to get the Butcher tableau for the third-order Heun's method, we have,
def BT(Method):
print(rk.loadRKM()[Method].dj_reduce())
BT('Heun33')
Heun RK 33
Heun's 3-stage, 3rd order
0 |
1/3 | 1/3
2/3 | 2/3
_____|_______________
| 1/4 0 3/4
Furthermore, we can prepare a Python code using the above algorithm.
def Heun_Method_3rd(f, y0, a, b, h= False, N=False):
'''
Inputs:
f: the ODE y'=f(t,y(x))
y0: the initial value
x_range: interval [a,b]
N: Number of points
h: step size h
'''
if N:
h = (b-a)/(N)
if h:
N = int((b-a)/h)
t = np.linspace(a, b, N+1)
y = np.zeros(t.shape, dtype=float)
y[0] = y0
for i in range(N):
k1 = f(t[i], y[i])
k2 = f(t[i]+ h/3, y[i]+ h*k1/3)
k3 = f(t[i]+ 2*h/3, y[i] + h*2*k2/3)
y[i+1] = y[i] + h*(k1 + 3*k3)/4
del k1, k2, k3
Table = pd.DataFrame({'t':t, 'y':y})
return Table
Example: We can use Heun's third-order method on the previous example IVP.
Header(Text = """Numerical Solution using the third-order Heun's method and N = 10""", C1 = Back.GREEN, C2 = Fore.GREEN)
Figs1(Solver = Heun_Method_3rd, Title = """Third-order Heun's method""")
Header(Text = 'Order of Accuracy')
Figs2(Solver = Heun_Method_3rd, h = [2**(-i) for i in range(4, 13)], yLim = [2.9, 3.1],
Title = """Third-order Heun's method""")
Line()
Numerical Solution using the third-order Heun's method and N = 10 ======================================================
| t | y | Exact | Error | |
|---|---|---|---|---|
| 1 | 0.1000 | 0.9950 | 0.9950 | 8.9306e-09 |
| 2 | 0.2000 | 0.9800 | 0.9800 | 4.1287e-08 |
| 3 | 0.3000 | 0.9551 | 0.9551 | 2.5240e-08 |
| 4 | 0.4000 | 0.9203 | 0.9203 | 2.7315e-07 |
| 5 | 0.5000 | 0.8762 | 0.8762 | 1.2114e-06 |
| 6 | 0.6000 | 0.8233 | 0.8233 | 3.1339e-06 |
| 7 | 0.7000 | 0.7627 | 0.7627 | 6.1831e-06 |
| 8 | 0.8000 | 0.6960 | 0.6960 | 1.0131e-05 |
| 9 | 0.9000 | 0.6250 | 0.6250 | 1.4309e-05 |
| 10 | 1.0000 | 0.5518 | 0.5518 | 1.7677e-05 |
Order of Accuracy ======================================================================================================
========================================================================================================================
This method is defined as follows, \begin{align} y_{n+1} &= y_{n}+\frac{h}{9}\left( 2k_{1} + 3k_{2} + 4k_{3} \right) \end{align} with \begin{align} \begin{cases} k_{1} & = f\left(t_{n},~y_{n} \right),\\ k_{2} & = f\left(t_{n}+\frac{h}{2},~y_{n}+\frac{h}{2}k_{1} \right),\\ k_{3} & = f\left(t_{n}+\frac{3}{4}h,~y_{0} + \frac{3h}{4}k_{2}\right) \end{cases} \end{align} or using the Butcher tableau \begin{align} \begin{array}{c|ccc}0&0&0&0\\1/2&1/2&0&0\\3/4&0&3/4&0\\\hline &2/9&1/3&4/9\\\end{array} \end{align}
Furthermore, we can prepare a Python code using the above algorithm.
def Ralston_Method_3rd(f, y0, a, b, h= False, N=False):
'''
Inputs:
f: the ODE y'=f(t,y(x))
y0: the initial value
x_range: interval [a,b]
N: Number of points
h: step size h
'''
if N:
h = (b-a)/(N)
if h:
N = int((b-a)/h)
t = np.linspace(a, b, N+1)
y = np.zeros(t.shape, dtype=float)
y[0] = y0
for i in range(N):
k1 = f(t[i], y[i])
k2 = f(t[i]+ h/2, y[i]+ h*k1/2)
k3 = f(t[i]+ 3*h/4, y[i] + h*3*k2/4)
y[i+1] = y[i] + h*(2*k1 + 3*k2 + 4*k3)/9
del k1, k2, k3
Table = pd.DataFrame({'t':t, 'y':y})
return Table
Example: We can use Ralston's third-order method on the previous example IVP.
Header(Text = """Numerical Solution using the third-order Ralston's method and N = 10""", C1 = Back.GREEN, C2 = Fore.GREEN)
Figs1(Solver = Ralston_Method_3rd, Title = """Third-order Ralston's method""")
Header(Text = 'Order of Accuracy')
Figs2(Solver = Ralston_Method_3rd, h = [2**(-i) for i in range(4, 13)], yLim = [2.9, 3.1],
Title = """Third-order Ralston's method""")
Line()
Numerical Solution using the third-order Ralston's method and N = 10 ===================================================
| t | y | Exact | Error | |
|---|---|---|---|---|
| 1 | 0.1000 | 0.9950 | 0.9950 | 2.1207e-06 |
| 2 | 0.2000 | 0.9800 | 0.9800 | 4.3623e-06 |
| 3 | 0.3000 | 0.9551 | 0.9551 | 6.7984e-06 |
| 4 | 0.4000 | 0.9203 | 0.9203 | 9.5125e-06 |
| 5 | 0.5000 | 0.8762 | 0.8762 | 1.2538e-05 |
| 6 | 0.6000 | 0.8233 | 0.8233 | 1.5762e-05 |
| 7 | 0.7000 | 0.7627 | 0.7627 | 1.8830e-05 |
| 8 | 0.8000 | 0.6960 | 0.6960 | 2.1097e-05 |
| 9 | 0.9000 | 0.6250 | 0.6250 | 2.1679e-05 |
| 10 | 1.0000 | 0.5518 | 0.5518 | 1.9596e-05 |
Order of Accuracy ======================================================================================================
========================================================================================================================
This method is defined as follows, \begin{align} y_{n+1} &= y_{n}+\frac{h}{9}\left( 2k_{1} + 3k_{2} + 4k_{3} \right) \end{align} with \begin{align} \begin{cases} k_{1} & = f\left(t_{n},~y_{n} \right),\\ k_{2} & = f\left(t_{n}+h,~y_{n}+hk_{1} \right),\\ k_{3} & = f\left(t_{n}+\frac{h}{2},~y_{0} + \frac{h}{4}k_{1} + \frac{h}{4}k_{2}\right) \end{cases} \end{align} or using the Butcher tableau \begin{align} \begin{array}{c|ccc}0&0&0&0\\1&1&0&0\\1/2&1/4&1/4&0\\\hline &1/6&1/6&2/3\\\end{array} \end{align}
The Butcher tableau can be also generated using nodepy python package. For example, to get the Butcher tableau for the third-order Heun's method, we have,
BT('SSP33')
SSPRK 33
The optimal 3-stage, 3rd order SSP Runge-Kutta method
0 |
1 | 1
1/2 | 1/4 1/4
_____|_______________
| 1/6 1/6 2/3
| 0.291 0.291 0.417
Furthermore, we can prepare a Python code using the above algorithm.
def SSPRK3(f, y0, a, b, h= False, N=False):
'''
Inputs:
f: the ODE y'=f(t,y(x))
y0: the initial value
x_range: interval [a,b]
N: Number of points
h: step size h
'''
if N:
h = (b-a)/(N)
if h:
N = int((b-a)/h)
t = np.linspace(a, b, N+1)
y = np.zeros(t.shape, dtype=float)
y[0] = y0
for i in range(N):
k1 = f(t[i], y[i])
k2 = f(t[i]+ h, y[i]+ h*k1)
k3 = f(t[i]+ h/2, y[i] + h*k1/4 + h*k2/4)
y[i+1] = y[i] + h*(k1 + k2 + 4*k3)/6
del k1, k2, k3
Table = pd.DataFrame({'t':t, 'y':y})
return Table
Example: We can use the SSPRK3 method on the previous example IVP.
Header(Text = """Numerical Solution using the SSPRK3 and N = 10""", C1 = Back.GREEN, C2 = Fore.GREEN)
Figs1(Solver = SSPRK3, Title = """SSPRK3""")
Header(Text = 'Order of Accuracy')
Figs2(Solver = SSPRK3, h = [2**(-i) for i in range(4, 13)], yLim = [2.9, 3.1],
Title = """SSPRK3""")
Line()
Numerical Solution using the SSPRK3 and N = 10 =========================================================================
| t | y | Exact | Error | |
|---|---|---|---|---|
| 1 | 0.1000 | 0.9950 | 0.9950 | 8.1570e-06 |
| 2 | 0.2000 | 0.9800 | 0.9800 | 1.5254e-05 |
| 3 | 0.3000 | 0.9550 | 0.9551 | 2.0552e-05 |
| 4 | 0.4000 | 0.9203 | 0.9203 | 2.3713e-05 |
| 5 | 0.5000 | 0.8761 | 0.8762 | 2.4931e-05 |
| 6 | 0.6000 | 0.8232 | 0.8233 | 2.4954e-05 |
| 7 | 0.7000 | 0.7627 | 0.7627 | 2.5002e-05 |
| 8 | 0.8000 | 0.6960 | 0.6960 | 2.6573e-05 |
| 9 | 0.9000 | 0.6250 | 0.6250 | 3.1160e-05 |
| 10 | 1.0000 | 0.5518 | 0.5518 | 3.9928e-05 |
Order of Accuracy ======================================================================================================
========================================================================================================================